#include "StdAfx.h"
#include "compute_matrix.h"
#include "matrix.hpp"
#include "evector.hpp"
#include "superelement.hpp"
#include "functions.hpp"
#include "Intergration.hpp"

#define LeftPoint 0
#define RightPoint 10
compute_matrix::compute_matrix(int ip_n)
{
	N = ip_n;
}

void compute_matrix::getStiffnessMatrix(FullMatrix<double> op_mt, vector<FElement> ip_vt)
{
    double **B;
	op_mt = FullMatrix<double>(N, N, 0);
    B = new double*[N];
    for(int i = 0; i < N; i++)
    {
        B[i] = new double[N];
    }

    for(int i = 0; i < N; i++)
        for(int j = 0; j < N; j++)
            B[i][j] = 0;

    for(int i = 0; i < N; i++)
    {
		int i0 = ip_vt[i].vertex0_index_;
		int i1 = ip_vt[i].vertex1_index_;
        B[i0][i0] += ip_vt[i].Bkk();
        if (i1<N)
		{
			B[i1][i1] += ip_vt[i].Bk1k1();
			B[i1][i0] += ip_vt[i].Bkk1();
		}
        B[i0][i1] += ip_vt[i].Bkk1();
        if(i0 == LeftPoint)
            B[i0][i0] += sigma0;
        if(i1 == RightPoint)
            B[i1][i1] += sigma1;
    }

    for(int i = 0; i < N; i++)
        for(int j = 0; j < N; j++)
            op_mt.Set(i, j, B[i][j]);
	op_mt.print();
}
void compute_matrix::getFVector(double *op_vt, vector<FElement> ip_vt)
{
    double *F;
	double d[2];
    F = new double[N];


    for(int i = 1; i < N; i++)
    {
		op_vt[i] = 0;
		int i0 = ip_vt[i].vertex0_index_;
		Intergration_Fk1(ip_vt[i], d, N);
        F[i0] += d[0];

		i0 = ip_vt[i].vertex1_index_;
		Intergration_Fk(ip_vt[i], d, N);
        F[i0] += d[0];
    }
	op_vt = F;
	for(int i = 0; i < N; i++)
	{
		std::cout << op_vt[i] << "\n";
	}
	printf ("\n");
}

compute_matrix::~compute_matrix(void)
{
}
